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Abstract 

Raman spectroscopy is the prime non-destructive characterization tool for graphene and 
related layered materials. The shear (C) and layer breathing modes (LBMs) are due to rela¬ 
tive motions of the planes, either perpendicular or parallel to their normal. This allows one 
to directly probe the interlayer interactions in multilayer samples. Graphene and other two- 
dimensional (2d) crystals can be combined to form various hybrids and heterostructures, cre¬ 
ating materials on demand with properties determined by the interlayer interaction. This is 
the case even for a single material, where multilayer stacks with different relative orientations 
have different optical and electronic properties. In twisted multilayer graphene samples there 
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is a significant enhancement of the C modes due to resonance with new optically allowed elec¬ 
tronic transitions, determined by the relative orientation of the layers. Here we show that this 
applies also to the LBMs, that can be now directly measured at room temperature. We find that 
twisting has a small effect on LBMs, quite different from the case of the C modes. This implies 
that the periodicity mismatch between two twisted layers mostly affects shear interactions. Our 
work shows that Raman spectroscopy is an ideal tool to uncover the interface coupling of 2d 
hybrids and heterostructures. 

Keywords: twisted multilayer graphene, layer breathing modes, interface coupling,first- 
principles calculations, resonant Raman spectroscopy,two-dimensional materials, two-dimensional 
heterostructures. 
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Layered materials ean be assembled to form heterostruetures held together by van der Waals 
interaetions. For a given assembly, the relative orientation of the individual layers ean ehange 
the optieal and eleetronie properties. This is also the ease when a single material is eonsid- 
ered. In multilayer graphene (MLG) samples, for a given number of layers (N), a wide range 
of properties is aeeessible by ehanging the relative orientation of the individual layers, 
refer to these as twisted-MLG (tMLG),^^ to indieate a mutual orientation of the planes differ- 
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ent from the naturally occurring one, with a twist angle The twist vector {p,q) is de¬ 

fined as the lattice vector of a supercell having q, p coordinates with respect to the basis vec¬ 
tors of single layer graphene (SLG). The twist angle can be derived from the twist vector as: 
cos0t={q^+4qp+p^)l2{q^+qp+p^). 

By assembling Bernal stacked m-layer (mLG, m > 1) and n-layer (nLG, n > 1) flakes, a 
(m-f-n)-system is formed, which we indicate as t(m-t-n)LG. In this notation, a Bernal-stacked 
BLG is denoted as 2LG, while a twisted one as t(l-i-l)LG. A flake consisting of a Bemal-stacked 
BLG placed at a generic angle Ot on a Bernal-stacked three layer graphene (TLG) is indicated as 
t(2-i-3)LG. This has significantly different properties when compared to a Bernal-stacked 5LG, or 
to a t(l-i-4)LG, or t(l-i-l-i-3)LG, etc, even though all these have the same N=5. For a given total N, 
the choice of m,n, etc. (with m-i-n-i-.. =N) and relative angles between each m,n,...LGs leads to a 
family of systems with different optical and electronic properties. Probing the coupling between 
the interface layers of mLG and nLG in t{m -\- n)LGs, and its impact on band structure and lattice 
dynamics, is crucial to gaining fundamental understanding of these systems and to tuning them for 
novel applications. 

Raman spectroscopy is one of the most used characterization techniques in carbon science 
and technology. The Raman spectrum of graphite and MLG consists of two fundamentally dif¬ 
ferent sets of peaks. Those, such as D, G, 2D, etc, present also in SLG, and due to in-plane 
vibrations, and others, such as the shear (C) modesand the layer breathing (LB) modes 
(LBMs),^®’^^’^^ due to relative motions of the planes themselves, either perpendicular or parallel 
to their normal. In NLG, all vibrational modes split due to the confinement in the direction per¬ 
pendicular to the basal plane, z, and, for a given N, there are N-1 C or LB modes, which we denote 
as Cnn-i and LBMa^a^-,- (/ = 1,2,...,A — 1), respectively. Here, Cni and LBMAfi(z.e., i = N —\) 
are the C and LB modes with the highest frequencies, respectively. However, due to the low elec¬ 
tron phonon coupling (EPC) and different symmetry, it has been not possible, thus far, to detect 
LBMs for samples at room temperature, unlike the highest energy C modes that can be measured 
in Bemal-stacked samples at room temperature. In Ref. 11 we have shown that, by performing 
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multi-wavelength Raman spectroseopy on tMLGs, an energy window exists, where a signifieant 
intensity enhaneement of the C peaks happens, due to resonance with new optically allowed elec¬ 
tronic transitions, determined by the relative orientation of the layers. This resonance effect is 
confirmed by the twist-angle dependence of the G and 2D intensities. 

Here we directly measure the LBM in tMLGs at room temperature with multi-wavelength Ra¬ 
man spectroscopy, and confirm their assignment by symmetry and polarization analysis combined 
with density functional theory (DFT). Similar to the C modes, the LBMs exhibit a significant 
intensity enhancement determined by the relative orientation of the layers. However, unlike the 
C modes, the observed LBMs are mainly determined by N, which suggests that the breathing 
coupling at the tMLG interfaces is almost independent of the relative layer orientation. The exper¬ 
imental positions of all LBMs can be described by a linear chain model considering next-nearest 
interlayer interactions, as verified by DFT. A charge density analysis reveals that the different be¬ 
havior of C and LB modes in tMLGs is due to the in-plane periodicity mismatch at the twisted 
interface. 


Results and Discussion 

The twisted samples are prepared as follows. Highly oriented pyrolytic graphite (HOPG) is me¬ 
chanically exfoliated on a Si/Si02 substrate. During exfoliation mLG flakes are folded onto 
nLG flakes to form t(m-i-n)LG flakes, such as the t(l-i-l-i-l)LG, t(l-i-3)LG, t(3-i-3)LG, t(4-i-4)LG and 
t(5-i-5)LG used in this study. Alternatively, a mLG flake from one substrate can also be trans¬ 
ferred onto a nLG flake on another substrate to form t(m-i-n)LG. Samples t(l-i-2)LG, t(2-i-2)LG and 
t(2-i-3)LG are prepared in this way. We follow the transfer method described in Ref.^^ A flake is 
exfoliated onto a polymer stack consisting of a water-soluble layer (Mitsubishi Rayon aquaSAVE) 
and PMMA, and the substrate is floated on the surface of a deionized water bath. During transfer, 
the target substrate is heated to 110°C to drive off any water adsorbed on the sample surface, as 
well as to promote good adhesion of PMMA to the target substrate. N in all initial and twisted 
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Figure 1: (a) Stokes/anti-Stokes Raman spectra in the C and LB spectral range, and Stokes Raman 
spectra in the G peak region for 1.96 and 2.33eV excitation. Polarized Raman spectra are also 
shown, (b) Peak area of C 31 , C 32 ,LBM 4 i and LBM 42 as a function of excitation energy. Solid 
diamonds, open squares and triangles are the experimental data, and solid and dashed lines are the 
simulations. The peak area of the Ei mode at 127 cm^^ of quartz, Aq^{E\), is used to normali z e 
all peaks. (c)A(LBM 4 i) as a function of excitation polarization direction. Open triangles are 
experimental data and solid lines are the expected trends the symmetry analysis. 
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MLGs is identified by Raman speetroscopy and optieal eontrast. 

Raman spectra are measured in back-scattering at room temperature with a Jobin-Yvon HR800 
Raman system, equipped with a liquid-nitrogen-cooled charge-coupled device (CCD), a lOOx 
objective lens (NA=0.90) and several gratings. The excitation energies are 1.58 and l.VleV from a 
Ti:Saphire laser, 1.96, 2.03, 2.09 and 2.28eV from a He-Ne laser, 1.83, 1.92, 2.18,2.34 and 2.41eV 
from a Kr+ laser, and 2.54, 2.67eV from an Ar+ laser. A 1800 lines/mm grating enables us to have 
each pixel of the charge-coupled detector cover 0.54cm^^ at 488nm. Plasma lines are removed 
from the laser signals using BragGrate Bandpass filters, as described in Ref. 21. Measurements 
down to 5cm~^ for each excitation are enabled by three BragGrate notch filters with optical density 
3 and with full width at half maximum (FWHM)=5- lOcm^ ^The typical laser power is~0.5mW 
to avoid sample heating. The accumulation time for each spectrum is~600s. 

We first consider a t(l-i-3)LG measured at 1.96 and 2.33eV, as for Figure 1(a). This shows 
peaks at~1510 and~1618cm“^ We assign these to the R and R' modes as described in Refs. 
30,31. From their position we deduce a 0/ ~10.6° between the SLG and TLG in this t(l-i-3)LG, 
see Methods for details. This corresponds to a twist vector (1,9). Two C modes (C 31 and C 32 ) 
are observed in t(l-i-3)LG, mainly localized in 3LG constituent, as previously discussed. Two 
additional modes are observed in t(l-i-3)LG at~l 16 and ~93cm~^ 

For a given N, the LBM position, Pos(LBM)n, can be written as:^°’^^ 


Pos(LBM)NN.i 


Pos(LBM)co sin 


in 

m 


( 1 ) 


where Pos(LBM)oo is the LBM in bulk graphite~128cm^^^^ We note that the N-1 LBM fre¬ 
quencies predicted by Eq. (1) do not necessarily translate to the experimental observation of the 
corresponding C and LBM Raman peaks, as these become Raman active under specific selection 
rules and symmetry constraints, as discussed in Methods. 

From Eq. (1) we get Pos(LBM2i)=90.5cm~^ and Pos(LBM 3 i)=l 10.8cm^^ The experimental 
value 116cm^^ is, however, larger than the predicted Pos(LBM 3 i), but closer to Pos(LBM 4 i)=l 18cm^^. 
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This implies that the LBM is eonsistent with that of a 4LG, but not with that of the 3LG eonstituent 
in the t(l+3)LG, unlike the case of the C modes, where the observed peaks correspond to C 31 and 
€ 32 ,^^ as indicated in 1. Thus, we assign the two LBMs in t(l+3)LG as LBM 41 and LBM 42 . 
Unlike the D and 2D modes, the LBMs are non-dispersive with excitation energy, Eex, as shown 
in Figure 1(a). This is expected, since they come from the Brillouin zone (BZ) center. The peak 
area of LBM 41 , A(LBM 4 i) measured at 1.96eV is~30 times higher than at 2.33eV, indicative of 
a resonant Raman behavior. We assign the LBM 41 and LBM 42 enhancement to resonance with 
new optically allowed electronic transitions in t(l+3)LG, as in the case of the C and G modes dis¬ 
cussed in Ref. 11. The C and LB modes are normalized to the Ei mode of quartz. Its position 
(~127cm~^) is so small that the CCD efficiency difference between C, LB and Ei modes for each 
excitation energy can be ignored. The resonant profile of EBM 41 is almost identical to that of C 32 , 
and the profile of EBM 42 is similar to that of C 31 , as shown in Eigure 1(b). This indicates that the 
EBM 41 resonant behavior can be also assigned to the resonance between the van Hove singularities 
in the joint density of states of all optically allowed transitions in t(l-i-3)EG and the laser excitation 
energy, similar to the C modes in tMEGs. 

Eigure 1(a) shows that the C and G modes are present in both parallel (XX) and cross (XY) 
polarization. However, the EBMs in t(l-i-3)EG vanish in the XY configuration. This can be ex¬ 
plained as follows. A t(m-l-n)EG (m 7 ^ n) has a C 3 symmetry, and the corresponding irreducible 
representation^^ is Y=A+E. All EBMs have A symmetry, all of C modes have E symmetry, and 
both the A and E modes are Raman active. The A Raman tensor is:^^ 


A = 


a 0 0 
0 a 0 
Q Q b 


( 2 ) 


This implies that, in backscattering, all EBMs should not be seen in the XY configuration, see 
Methods, and that their intensity is a function of the angle ((^) between the polarization of the 
incident light and the polarization (Y) of the Raman signal, I (LBM) = a^cos{(j))^ (see Methods). 
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Figure 1(c) plots /(LBM 41 ) as a function of 0. The experimental data (open triangles) are in good 
agreement with the symmetry analysis. 



Figure 2: Stokes/anti-Stokes Raman spectra in the C and LB peak region and Stokes spectra in the 
G spectral region for six tMLGs. ^ex is also indicated. The spectra are scaled and offset for clarity. 
The scaling factors of the individual spectra are shown. Vertical lines are guides to the eye. 

Figure 2 plots the Raman spectra of six tMLGs: t(l+2)LG, t(l+l+l)LG, t(l+3)LG, t(2+2)LG, 
t(2+3)LG and t(5+5)LG. To facilitate comparison, all are normalized to have the same intensity 
of the G peak, 1(G). The spectra show the C modes of mLG (m >1) and nLG (n >1), localized 
inside the mLG or nLG constituents.^^ However, this it is not the case for the LBMs. E.g., in 
t(l+l+l)LG there is no observable C mode, because the twisted interface significantly weakens 
the shear coupling and pushes the C frequency towards the Rayleigh line, outside the measured 
spectral region. However, in the LBM region, t(l+l+l)LG shows a peak at~108.8cm^\ close to 
the predicted LBM 31 ~110.8cm^^ A similar peak at ~109.9cm^^ is observed in t(l+2)LG. Since 
both t(l+l+l)LG and t(l+2)LG are two possible t3LG embodiments, we assign the two LBMs 
in t(l+l+l)LG and t(l+2)LG to LBM 31 . The t(2+2)LG sample shows a LBM~115.5cm^\ very 
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close to the observed~116em ^ in t(l+3)LG, and to the expeeted value for LBM 41 . However, 
unlike t(l+3)LG, t(2+2)LG has a D 3 symmetry, and LBM 41 and LBM 43 are Raman-aetive Ai 
modes, while LBM 42 is a Raman-inaetive A 2 mode, see Methods. Thus, LBM 42 in t(2+2)LG is 
not deteeted due to symmetry. In a similar way, we assign the LBMs in t(2+3)LG and t(5+5)LG as 
LBM 51 , LBM 52 and LBMioj, respeetively. Based on symmetry, all C modes in t(m + n)LGs are 
Raman aetive. Consequently, the C modes of the Bernal-staeked eonstituents are also observed, 
sueh as C 51 , C 53 and C 54 in t(5+5)LG. 

The above data suggest that, unlike the the C modes, Pos(LBMAf^iv-!) in a tAfLG {N = m-\-n-\- 
...) is mainly determined by N and not by the number of layers of the individual Bernal-staeked 
eonstituents (m,n,...). This means that the LBMs in tMLG are not loealized inside its eonstituents, 
but are a eolleetive motion involving all layers. We stress that dt for the six tMLGs in 2 is not the 
same, as determined by the respeetive R' and R positions. Various dt give different band struetures 
with different values for optieally-allowed resonance transitions. Therefore, for each sample 
we deteet LBMs at different exeitations. 

We now eonsider the effeets of ehanging interlayer interaetions on the LBM positions. To do 
so, we solve the equation of motion for a linear ehain system.The frequeneies G) (in em^^) and 
displaeement patterns ean be ealeulated by solving linear homogeneous equations: 


(ofm = 


^ ^ ^ Du/, 


(3) 


where u,- is the phonon eigenveetor of the ith mode with frequeney ( 0 /, /i=7.6x 10“^^kgA“^ is the 
SLG mass per unit area, c=3.0x lO^^em s~^ is the speed of light, and D is the foree eonstant matrix. 
In our previous works, we adopted a simple linear ehain model (LCM) with only nearest-neighbor 
interlayer interaetions. This allowed us to explain the observed C modes in Bernal and tMLGs, 
as well as the LBMs in several 2d materials.For tMLGs, this also prediets the C modes 
by introdueing a weaker shear foree eonstant (aj^) at the twisted interfaee. 

The top panel of Figure 3(a) plots the sehematie LCM for LBMs in t(2-i-3)LG if only the 
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Figure 3: (a)Linear chain model (LCM) and LCM with second-nearest interlayer coupling 
(2LCM). (b) Theoretical (LCM, open triangles; 2LCM, open diamonds) Pos(LBMAri) and 
Pos(LBMAf 2 ) in 4LG and 5LG, and experimental (Exp., crosses) and theoretical (2LCM, open 
diamonds) Pos(LBM 7 vi) and Pos(LBMa? 2 ) in t(2+2)LG, t(l+3)LG and t(2+3)LG. (c) Experimental 
(Exp., open crosses) and theoretical (2ECM, open diamonds) Pos(LBMjvi) and Pos(LBMiV 2 ) in 
tNLG. (d) Normal mode displacements and frequencies of t(l+3)EG and t(2+3)LG based on the 
2LCM. 
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nearest-neighbor layer-breathing interlayer interaction is considered. The experimental fre¬ 
quencies of t(2-i-2)LG, t(l-i-3)LG and t(2-i-3)LG are plotted in Figure 3(b) as crosses, and those of 
all tMLGs are summarized in Figure 3(c), including LBMs from t(3-i-3)LG and t(4-i-4)LG, whose 
Raman spectra are presented in Methods. By taking the average frequency (115.8cm^^) of the 
experimental LBM 41 measured in t(l-i-3)LG and t(2-i-2)LG, we get Ci5q“=106x 10^® Nm~^, which 
would give 119.2cm^^ for Pos(LBM 5 i) in 5LG, consistent with the value measured in t(2-i-3)LG. 
Figure 3(b) also gives Pos(LBM 42 )= 88 . 6 cm^^ for 4LG and Pos(LBM52)=101.4cm~^ for 5LG, 
which are 4.3 and 2.9cm~^ lower than those observed in t(l-i-3)LG and t(2-i-3)LG, respectively. 

These lower frequencies suggest that the LCM, with only nearest-neighbor interlayer interac¬ 
tions, may be insufficient to reproduce the interlayer breathing coupling in tMLGs. If a weakened 
coupling at the twisted interface is included in the LCM, it will result in LBM red-shift for both 
LBMjvi and LBMa ?2 (N=4,5), making the agreement worse, see Methods. We thus introduce an 
interlayer breathing force constant between the second-nearest neighbor layers (/3 q^). The new 
model is denoted as 2LCM, and is schematically shown in Figure 3(a) for LBMs in t(2-i-3)LG. 
2LCM with /3 ,^~9.3xl0^^ Nm ^ fits the experimental data best, as indicated by diamonds in 
Figure 3(b). With 2LCM we can well fit the frequencies of the observed LBMs in all tMLGs, as 
shown in Figure 3(b,c). Additionally, we can expand the 2LCM predictions to bulk graphite, based 
on the parameters fitted on our experiments. The silent LBM (B 2 g) in graphite is derived to be 
~125.3cm^^, very close to ~128cm~^ determined by neutron spectrometry.^^ 

The normal mode displacements and frequencies of each LBM in t(l-i-3)LG and t(2-i-3)LG as 
derived by the 2LCM are summarized in Figure 3(d). In LBM/v 1 , the relative motions of the 
nearest-neighbor layers are always out-of-phase, and those of the second-nearest-neighbor layers 
are always in-phase. This would suggest Pos(LBMiv,i) to be insensitive to the second-nearest- 
neighbor interlayer coupling. However, the relative motions of the second-nearest-neighbor are 
out-of-phase for LBM42 in t(l-i-3)LG and LBM52 in t(2-i-3)LG. Thus, the reason why Eq. ( 1 ) fits 
Pos(LBM 7 vi) well, but predicts lower frequencies for Pos(LBMiV2) is, most likely, due to the lack 
of interaction from second-nearest-neighbor layers. 
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The 2LCM gives the same LB eoupling for twisted and Bernal-staeked interfaees. However, 
the shear eoupling at twisted interfaees is~20% of that at Bernal-staeked interfaees. We now use 
DFT and density funetional perturbation theory (DFPT)^^ to validate this model, and to understand 
the differenee between the C and LBMs in tMLG. Beeause a t(m-t-n)LG with a twist veetor of 
(1,2), i.e., a twist angle of 21.8°, is a simplest twist strueture, we eonsider t(2-i-3)LG and t(l-i-2)LG 
with this twist angle for DFPT. 

Table 1: Ab initio interlayer foree eonstants between eaeh eouple of layers along z for t(2-i-3)LG. 
Twisting happens between the seeond (denoted 2a) and third (denoted 3b) layers. Two eategories 
of Bernal-staeked layers are grouped as "a" and "b", respeetively. 


Foree eonstant 
(xlQi^Nm-^) 

la 

2 a 

3b 

4b 5b 

la 

- 

- 

- 

- 

2 a 

114.2 

- 

- 

- 

3b 

4.3 

114.7 

- 

- 

4b 

3.9 

3.4 

120.1 

- 

5b 

4.0 

6.1 

2.6 

113.4 


We first ealeulate the frequeneies of LBMs in t(2-i-3)LG with a (1,2) twist veetor. They are 
126.3em~i (LBM 51 ), 107.3em"i (LBM 52 ), 19.9cm ^ (LBM 53 ) and 47.9em-i (LBM 54 ), respee¬ 
tively, overall eonsistent, but a few em^^ larger, than the experiments reported in Figure 3, owing 
to the slightly overestimated interlayer interaetion. A full eomparison between ealeulated and 
measured frequeneies is reported in Methods. Figure 3(a) shows the five layers and four inter¬ 
faees in t(2-i-3)LG. We denote them as la, 2a, 3b, 4b, and 5b from left to right. Twisting happens 
between layers 2a and 3b, and we eall this interfaee 2a-3b. The interlayer foree eonstant (IFC) 
along z is a measure of the interlayer breathing eoupling and is ealeulated as for Methods. The 
IFC along z between 2a and 3b (the twisted interfaee) is 114.7 x lO^^Nm'^, elose to that of other 
Bernal staeked interfaees, pointing to a similar breathing eoupling at the twisted interfaee as that 
of the Bernal staeked interfaee. We get a(^=115.6x lO^^Nm^^ by averaging the eomputed IFCs 
along z at the four interfaees, whieh agrees with the experiments and with the value derived using 
the 2LCM ( 116 x lO^^Nm^^) 
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Figure 4: (a) Atomic structure of t(l+2)LG with 0^=21.8°. The twisted and Bernal-stacked inter¬ 
faces are indicated by green solid and dashed lines, respectively, (b) Charge density contour for 
the t(l-i-2)LG sample in a. (c) Layer-averaged charge density along z. 

We now address the substantial force constant difference for the C and LBMs in twisted and 
Bemal-stacked layers. Van der Waals forces, specifically the dispersion force,'^® rule the interlayer 
interactions, and play a key role in the difference between C and LBMs in tMLG. Figure 4(a) 
plots the sideview of the fully-relaxed atomic structure of a t(l-i-2)LG with a (1,2) twist vector. We 
also consider t(l-i-2)LGs with twist angles of 13.2°, 38.2° and 46.8°. Four stacking configuration 
are considered for each angle. The average interlayer distance for every configuration is 3.39A, 
with a variation less than O.OIA. Ref.'^^ reported a similar result for twisted M 0 S 2 bilayers, with the 
calculated interlayer distances nearly identical in the 0° to 60° range. Our calculations are also con¬ 
sistent with the interlayer distance in t(l-i-l)LGs calculated in Ref.,^^ showing a larger interlayer 
distance at the twisted interface when compared to Bemal-stacked layers, and little correlation 
between interlayer distance and twist angle. The interlayer distance between the twisted interface 
and the Bemal-stacked interface in t(l-i-2)LG is~0.1A , much smaller than in MoS 2 /MoSe 2 het- 
erostmctures (~0.6A),^^ where the interface has~4% lattice mismatch. This is directly relevant 
for the out-of-plane breathing vibration along z, as represented by the LBM frequency. Eq. (4) 
and Eq. (5) in Methods indicate that the interaction strength has a positive correlation with charge 
density, nearly identical at the two interfaces of Eigure 4(b). A small difference is revealed by 
calculating the mean charge densities at the two interfaces. The interlayer breathing interaction at 
the twisted interface is very close to that of Bemal-stacked interfaces, again supporting the 2ECM. 
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Figure 5: Atomic structure of (a) Bernal-stacked 2LG on the top of t(l-i-2)LG and (b) the t2LG 
at the bottom of t( 1-1-2)LG. The corresponding charge density (c) at the Bernal-stacked interface 
in (a) and (d) at the twisted interface in (b). Schematic diagram for the charge distribution (e) 
at the Bemal-stacked interface in (a) and (f) at the twisted interface in (b). The latter shows the 
mismatched periodicity between the two layers. 


We turn to consider the C modes in t(l-i-2)LG with a (1,2) twist vector. Top views of the 
Bemal-stacked and twisted interfaces are shown in Figure 5(a,b), while their corresponding charge 
densities in the middle of two SLGs is shown in Figure 5(c,d). Both plots indicate that the Cg 
symmetry at the Bemal-stacked interface is broken at the twisted interface (Figure 5(b)), and the 
local density periodicity is also lifted (Figure 5(d)). Twisting forms a Moire pattern, resulting in 
a locally mismatched periodicity of the charge density variations. Figure 5(e,f) plots a schematic 
diagram illustrating the effect of periodicity mismatch on the C vibrations. In Bemal-stacked inter¬ 
faces the interatomic restoring forces are all along the positive direction for a small displacement. 
Figure 5(e). With the elimination of the local periodicity, a Moire pattern at the twisted interface 
makes the interatomic restoring forces negative or positive, as shown in Figure 5(f). Therefore, 
shear restoring forces are nearly canceled at the twisted interface, resulting in a much weaker 
shear coupling than in Bemal-stacked interfaces. Thus, the softening of the C modes is due to the 
periodicity mismatch at the twisted interface. 
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Conclusions 


We measured by resonant Raman spectroscopy the LBMs of tMLG, an archetype heterostructure. 
We showed that a second-nearest neighbor linear chain model explains all the measured spectra, 
as validated by ab-initio calculations. The interlayer shear coupling strength declines at twisted in¬ 
terfaces due to the periodicity mismatch between two twisted layers, while the interlayer breathing 
coupling remains nearly constant. Beyond tMLGs, the interlayer interaction of other heterostruc¬ 
tures can also be measured by Raman spectroscopy."^^ Unlike graphene, the interlayer coupling 
modes of other 2d layered materials, like transition metal chalcogenides (e.g. M 0 S 2 and 

WSe 2 ) and others, such as NbSe 2 ,‘^^ and Bi 2 Se 3 "^"^ and Bi 2 Te 3 ,"^"^ can be measured more easily, 
due to the stronger electron-phonon coupling. Therefore, the LBMs should be also measurable in 
heterostructures with clean interfaces, such as graphene/MoS 2 , grapheneAVS 2 , MoS 2 AVSe 2 , thus 
allowing one to probe the interlayer coupling of these two-component layered heterostructures and, 
possibly, even more complex structures. By studying both C and LB modes together, it should be 
possible to detect the detailed components, number of layers of each component, and the coupling 
amongst the components, a crucial step for both fundamental science and technology based on 
these materials. 
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Methods 


Calculations 

Structural relaxation and charge density calculations are performed using the DFT code Vienna 
ab-initio simulation package (VASP)^^ within the projector augmented wave method"^^’^^ and a 
plane-wave basis. The exchange-correlation potential is treated within the generalized gradient 
approximation. Van der Waals interactions are considered under the framework of the vdW-DF 
method'^® with the optB86b exchange functional. This exchange-correlation combination is more 
accurate in predicting lattice parameters in 2d materials, such as black phosphorous^^ and boron 
nitridethan other vdW-DF approaches, while it is known to slightly overestimate interlayer 
binding energy. In vdW-DF the description of the dispersion force requires the inclusion of 
the non-local correlation energy 

Ef = ^ j y drdr'n(r)d>(r,r')n(r') (4) 

o 4 

d>fr —y --- (51 

2m2(Oo(r)(OD(r') [cOD(r) + tOD(r')]^^ 

with n{r) the charge density, $ the correlation interaction kernel and d the distance between two 
SLGs. For d ^ oo, ^ oc n~^'^d~^, which means Ef oc n^'^d~^. The non-local correlation energy 
between two SLGs is determined by charge density and layer distance. 

A 29 X 29 X1 k-mesh is used to sample the BZ for Bemal-stacked supercells and an 11x11x1 
one for twisted supercells, due to the y/l larger lattice constant. The energy cutoff for the plane- 
wave basis is 400eV. All atoms are fully relaxed until the residual force per atom is smaller than 
O.OOleV-A ^ Vibrational frequencies are calculated using DFPT,^^ as implemented in VASP. 
In an interlayer vibrational mode, the whole layer can be treated as one rigid body. The IFC is 
constructed by summing inter-atomic force constants over all atoms from each of the two adjacent 
layers. The matrix of inter-atomic force constants, essentially the Hessian matrix of the Bom- 
Oppenheimer energy surface, is defined as the energetic response to a distortion of atomic geometry 
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in DFPT. 


37 

Figure 6 plots the LBM positions ealeulated from the LCM in Eq. (1) and by DFT for Bernal 
Staeked samples (with DFT data rigidly shifted by ~10em^^) implemented in the QuantumE- 
SPRESSO paekage.^^ The in-plane lattiee eonstant is set to 2.43 A and the interlayer distanee 
3.26 A to match the experimental ZO' frequency at the F point. A norm-conserving Martins- 
Troullier pseudopotential within the local density approximation (EDA) is used, and the plane 
waves were expanded up to a 80Ry cutoff. The BZ is sampled using a 12x12x4 Monkhorst-Pack 
mesh and Methfessel-Paxton smearing with 0.03Ry width is used for the electronic occupations 
close to the Eermi level. The dynamical matrices are computed on a 8 x 8 x 3 mesh. The modes are 
either Raman (R) or Infrared (IR) active. 
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Eigure 6: EBMs for Bernal stacked NEGs. Gray lines indicate the calculated ECM. 

Eigure 7 compares DEPT and experimental Pos(C) and Pos(EBM) in various tMEGs. 

Eigure 8 compares the experimental EBMs in tMEGs with those calculated with the ECM of 
Eq. (1) and those using a ECM with a weakened coupling at the twisted interface (tECM). A 10% 
weakened coupling red-shifts both EBM^yi and EBM^f 2 (N=4 and 5), resulting in a worse fit to the 
experimental data. 
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Figure 7: C and LBMs in tMLGs. (pink crosses) DFPT data, (blue circles) Experiments. 
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Relation between Of and the frequency of R and R' modes 

The observation of R and R' peaks in the Raman speetra of tBLG is due to the superlattiee mod¬ 
ulation aetivating phonons in the BZ interior. Of dietates the waveveetor for this modulation, 
with q{9t) the difference between the basic vectors of two SLGs in the BZ. The waveveetor q{0t) 
selects the phonons along the phonon dispersion that become Raman active. The relation between 
qTK{9t) and the 9t is given by: 


qrK{9t) = f 1 ~ ^7 — 2\/hsin9t — 6cos9t ^ , 


3a 




( 6 ) 


where a=2.46A is the SLG lattice constant. From Pos(R) and Pos(R'), qrK{9t) can be determined 
from the SLG phonon dispersion. Eq. ( 6 ) then gives 9t. For the assignment, we use phonon dis¬ 
persions calculated from DFT^^ corrected with GW (Green's function G of the screened Coulomb 
interaction W), which well reproduces the experimental FO-TO splitting. 

Figure 9 plots the optical image and Raman spectra of t(l-i-l-i-l)FG and t(l-i-l)FG. There are 
two couples of R and R' modes in t(l-i-l-i-l)FG due to twice folding a SFG. The Ri mode of 
t(l-i-l-i-l)FG is at 1529cm \ the same position as the R mode of t(l-i-l)FG. This means that the 
Ri and Ri' are from the bottom twisted bilayer of t(l-i-l-i-l)FG and that R 2 is from the top twisted 
bilayer of t(l-i-l-i-l)FG. 

Figure 10 plots the optical image and Raman spectra of t(3-i-3)FG and t(4-i-4)FG. 9t of t(3-i-3)FG 
and t(4-i-4)FG are 11.4° and 12.0°, respectively, determined by the respective R' modes. 


Symmetry and Raman activity of C and LBMs in t(m + «)LG (m ^ n) and 

t(« + «)LG {n > 2) 

t(m-t-n)FG (m 7^ n) have C3 symmetry, the corresponding irreducible representation is r=A-i-£', 
and both A and E modes are Raman active.In t(m-t-n)FG with (m 7 ^ n), all non-degenerate 
FBMs have A symmetry, and all of double-degenerate C modes belong to E symmetry. 

t(n -f n)FG {n > 2) have D 3 symmetry, and the corresponding irreducible representation is 
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Figure 9: Optical image and Raman spectra of t(l+l+l)LG and t(l+l)LG. 



Figure 10: Stokes/anti-Stokes Raman spectra in the C peak region and Stokes spectra in the LBM 
and G spectral regions for t(3+3)LG and t(4+4)LG. The twist angle and laser energy is marked for 
each sample. 
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r=Ai+A 2 +£'.^^ Ai and E modes are Raman aetive, while A 2 are Raman inactive.In t(2+2)LG, 
LBM 41 and LBM 43 have A\ symmetry, while LBM 42 has A 2 symmetry, and all the C modes are 


The Raman intensity is proportional to \ei-Rt • where e, and are the unit vectors describ¬ 
ing the polarizations of the incident and scattered light, and Rt is Raman tensor. In our work, 
the polarization of the incident light is at an angle (0) set by a A /2 wave plate (e/=[cos0 sin 0 0]), 
and the polarization of the scattered light is fixed along the horizontal (e^=[l 0 0]')- Therefore, the 
Raman tensors of for the LBMs in t(m -f n)LG (m ^ n) is: 


A = 


a 0 0 
0 a 0 
0 0b 


Thus, I(LBMs) in t(m + n)LG (m 7 ^ n) is: 




a 

0 

0 


1 


l{LBMs) oc 

[cos^ sin^ 0 ] 

0 

a 

0 


0 




0 

0 

b 


0 



/y ^ 

a cos 


(7) 


( 8 ) 


As discussed above, the LBMs in the t(l-i-3)LG are Raman active, except LBM 42 (Raman inactive). 
LBM 41 and LBM 43 in t(2-i-2)LG are Raman active. Both LBM 41 and LBM 42 are observed in 
t(l-i-3)LG, see Figure 2. However, only LBM 41 is observed in t(2-i-2)LG. The absence of LBM 43 
in t(l-i-3)LG and t(2-i-2)LG may result from a weaker EPC.^^ The Raman tensor of the A 2 mode in 
t(2-i-2)LG is the same as that of the A mode in t(l-i-3)LG,^^ thus the I(LBM 4 i) in t(2-i-2)LG is also 
laser-polarization dependent. 
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